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Abstract 

The behavior of non-local thermal-equilibrium (NLTE) plasmas plays a central role in many fields of modern-day 
physics, such as laser-produced plasmas, astrophysics, inertial or magnetic confinement fusion devices, or X-ray 
sources. The proper description of these media in stationary cases requires to solve linear systems of thousands or 
more rate equations. A possible simplification for this arduous numerical task may lie in some type of statistical 
average, such as configuration or superconfiguration average. However to assess the validity of this procedure and to 
handle cases where isolated lines play an important role, it may be important to deal with detailed levels systems. This 
involves matrices with sometimes billions of elements, which are rather sparse but still involve thousands of diagonals. 
We propose here a numerical algorithm based on the LU decomposition for such linear systems. This method turns 
out to be orders of magnitude faster than the traditional Gauss elimination. And at variance with alternate methods 
based on conjugate gradients or minimization, no convergence or accuracy issues have been faced. Some examples 
are discussed in connection with the krypton and tungsten cases discussed at the last NLTE meeting. Furthermore, to 
assess the validity of configuration average, several criteria are discussed. While a criterion based on detailed balance 
is relevant in cases not too far from LTE but insufficient otherwise, an alternate criterion based on the use of a fictive 
configuration temperature is proposed and successfully tested. It appears that detailed calculations are sometimes 
necessary, which supports the search for an eflncient solver as the one proposed here. 
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1. Introduction 

Plasmas in a non local thermodynamic equilibrium (NLTE) state play an important role in several domains of 
physics, such as in astrophysics, in magnetic and inertial confinement devices, in radiation-solid target interaction 
experiments at large laser facilities, or in XUV and X-rays sources! ll 2 ] . The proper calculation of important quantities 
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such as the opacity or the emission efficiency requires an accurate atomic physics description of the active medium. In 
time-independent regimes this impHes to solve large systems of kinetic equations describing the population transfers 
between each level. Several codes have been developed to deal with NLTE plasmas, among which ATOMIC jQ, 
SCROLL |4L MOST isl, ATOM3R |Ql, AVERROES |7], FLYCHK |8], CRET 
SCAALP 1 3l- The NLTE meetings provide an opportunity to benchmark them | 

Because of its inherent complexity, NLTE physics, even in its stationary form, lies on efficient linear algebra 
algorithms. In our previous approach to this subject ifni . we relied on the standard Gauss algorithm: it proves to be 
very stable and accurate for systems up to several thousand equations, but it becomes prohibitively slow above. A 
possible workaround to this issue is to resort to statistical averages, the simplest being the configuration average (CA), 
or the superconfiguration average. A discussion of these averaging procedure efficiency may be found in the literature 
(e.g., jlSll ). However, using this average one may neglect important properties, for instance when the various levels 
in a given configuration have very different excitation rates. 

This paper investigates two complementary directions to allow the treatment of complex NLTE plasmas. First, we 
discuss various possible matrix-inversion algorithms, this step usually being the bottleneck of the collisional-radiative 
codes (Sec. |2|. Then, in Sec. [3] we review several possibilities aimed at designing a criterion that might efficiently 
qualify the CA procedure. Sees. |4] and |5] illustrate this discussion by dealing with two cases considered at the last 
NLTE meeting ildl . Brief concluding remarks are given in Sec.|6] 



2. Comparison of sparse-matrix inversion techniques 

2.1. System properties 

In order to properly characterize plasma properties, including observables as simple as the average net ion charge, 
one needs to account for a large set of ionic states, including series of doubly or even multiply excited states. This 
implies that the rate matrices describing the radiative and coUisional transfers in such plasmas may have a huge 
number of elements. Because usually the only transition processes considered in collisional-radiative systems induce 
a net charge variation AZ* = 0, +1, these matrices have a tridiagonal-block structure as shown in Fig.[T] If a lot of 
Z* are considered together, and if the numbers of levels considered for each ion are similar, the matrix contains a 
large proportion of zeroes and is thus reasonably sparse. Even more, inside the tri-diagonal blocs, some transition 
rates may cancel (such as those between multiply and singly excited states). However the algorithms proposed here 
do not account for such zeroes inside the tri-diagonal blocs. A direct estimate in a krypton case illustrated by Fig.[T] 
revealed that 30 to 40% of such elements do cancel: for this 5586 x 5586 matrix, with 3.1 x 10^ elements, 1.2 x 10^ 
are in diagonal blocs, which contain themselves 3.8 x 10^ (32%) "accidentally" null elements. However, dealing 
properly with such "accidental" cancelation requires to use a cumbersome indexing scheme which may slow down 
computations, and such attempt was not made in the present work. Clearly this sparse character is stronger when 
many charge states are included and when each bloc has approximately the same size. 
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2.2. Gauss algorithm 

The first algorithm used to solve collisional-radiative systems is the standard Gauss-Jordan elimination. We could 
check that this method provides in all considered cases a fair numerical accuracy: this can be done for instance by 
cancelin g the unbalanced transition rates and by controlling the level-population departure from the Saha-Boltzmann 
solution I IT*]. Except when one pivot accidentally cancel — which was never observed in properly defined NLTE 
cases — the algorithm always provide a solution, through a constant number of operations, which scales as the cube 
of the number of equations. To sum up, the Gauss elimination is robust, accurate, but it is really slow and does not 
take benefit from the sparse character of the system. It will be used only when we consider small systems — as those 
in configuration average discussed below — or when we wish to cross-check other algorithms. 

2.3. Conjugate gradients methods 

Conjugate Gradients (CG) methods are iterative procedures for solving a linear system through a series of ana- 
lytically defined steps. The biconjugate gradient method generalize the procedure to non symmetric definite positive 
matrices. The convergence of the process may be linked to the use of a preconditioner as discussed, e.g., in 
and will not be considered here. In this paper our goal is to check how standard routines widely available perform, 
evaluating their general robustness regardless of the many variants that exist. One must notice that, at variance with 
ifioi where non-stationary systems with few hundreds of levels were included, one may deal here with several tens 
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of thousands of equations. Using CG methods — here we use linbcg from Numerical Recipes 
from the sparse character of the rate matrix since the user is supposed to provide an efficient way to perform the matrix 
by column vector product. 

It appears that the CG codes with no preconditioner may poorly work. As seen below, such iterative methods start 
on a user-defined solution, and the convergence may strongly depend on this step. Besides a convergence criterion 
must be chosen. It was found that the inline criterion included with linbcg routine is not satisfactory: The proposed 
error parameter may be still about 1 while convergence has indeed be reached. Therefore in this work we used the 
natural and physically sound criterion based on the average net charge Z* variation for the last two iterations. Namely, 
when for the first time 

I < Z* >„ - < Z* >„_! I < s and | < Z* >„ - < Z* >„_2 \ <e (1) 

with for instance e - 10"^ convergence is assumed to be reached. Two steps are considered since it appeared that a 
single check could produce an artificial convergence. 

As shown below (subsection 12.51 1. the CG methods when they converge are fast and accurate. However, the CG 
algorithm may lack of robustness. For instance, we considered the case where radiative and autoionization rates 
are canceled: including only the collisional excitation and ionization and the reverse processes, and comparing the 
solution with Saha-Boltzmann is a way to check the algorithm accuracy In a krypton plasma at a 500 eV 

temperature and 10'** cm""* electronic density, with 5 662 levels and 14 charge states accounted for, we did obtain 



populations different from Saha-Boltzmann by 2 x lO"'-' in average (6 x 10"'*^ at maximum), but this was obtained 
after a rather large number of iterations: 1966 with s - 10 ' in Eq. ([TJ, starting with zero populations. The inspection 
of < Z* > before convergence shows strong fluctuations from one iteration to the next. One may estimate that the 
condition e - 10 is too severe, which overestimate the convergence effort, but if one stops the process at, e.g., 150 
iterations as obtained in well behaved cases (see next paragraphs), the < Z* > fluctuations on the last iterations are 
above 10"^, which is unacceptably large. 

Moreover, as for every iterative process, the CG convergence may depend on the initial choice for the populations. 
To illustrate this when no preconditioner is used, we have plotted in Fig. |2]and Fig. |3]results of a NLTE analysis in 
krypton for an electronic temperature of Tg = 1000 eV and an electronic density of A^g = lO''* cm-3. Here the number 
of equations is N = 5586, large enough to check significantly the calculation speed, but still moderate to allow direct 
computation by the slow but stable Gauss elimination. The linbcg internal convergence criterion gives an "error" 
of 3.3 after 140 iterations and 19.6 after 150 iterations, while it is obvious from Fig. [3] that convergence does then 
occur One notices that initializing the process with zero populations give fast and accurate convergence: 150 steps if 
accuracy of 10 on < Z* > is sought as indicated by the lower part of Fig.[3j The otherwise known correct < Z* > 
is given by the horizontal line in the upper figure. 

Conversely, initialization on Saha-Boltzmann populations leads to an erratic behavior. After about 4000 steps one 
gets close to the correct solution but with a mediocre 10^^-10"^ accuracy (cf. lower Fig.|2|i. Then the solution begins 
to diverge and reaches after almost 10 000 iterations a "false convergence" at < Z* >^ 28.162, significantly off the 
(otherwise confirmed) Gauss solution of < Z* >= 27.660. Up to 25 000 iterations have been performed here. The 
< Z* > variation from one step to the next as seen on the lower part of Fig. |2] remains above 10"^: The fact that 
this figure is much greater than machine accuracy ^ 10"' indicates that this "convergence" is not satisfactory. This 
further indicates that the s parameter of criterion ^ must be chosen very carefully. 

2.4. LU decomposition 

Another possible choice for matrix inversion is based on the LU decomposition algorithm, which uses an upper- 
lower triangular matrix decomposition. The used routine is here dgbsvx from lapack linear algebra package, which 
assumes a band-diagonal structure for the matrix to be inverted. The balancing of rows and columns is automatically 
done inside the code. The matrix condition is checked and given as an output parameter; In the various cases tested, 
at or far from thermal equilibrium, it was found to be satisfactory. In addition to the storage of the original matrix 
with Kl and Ku lower and upper diagonals, which occupies (Ki + Kjj + l)xN, N being the number of equations, the 
proposed algorithm requires the storage of the lower and upper triangular matrices, i.e., (2Ki + Ku + \)y.N additional 
memory locations. Considering reaUstic values for the (- Ku) parameter, this amounts to foresee storage as large 
as 32 GB for a matrix size ^ 50 000. 

The various tests performed have shown that this algorithm is fast, accurate and robust. No particular row or col- 
umn balancing had to be done by a preprocessor, since this is automatic in the lapack routine. The matrix condition 
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is checked inside the code and it was found satisfactory in every considered cases. As an accuracy test proposed 
before ifni . we solved with the LU algorithm a kinetic system of 448 1 equations including only collisional excitation, 
collisional ionization plus the reverse processes. One gets then the Saha-Boltzmann solution with a maximum pop- 
ulation difference of 1.2 x 10"'^, and an average difference of 5.4 x 10"'*, the Gauss algorithm providing the same 
excellent accuracy. It appears that the only serious limitation to LU method is the storage issue, even for reasonably 
sparse matrices. Looking again at Fig.[Tl the stored elements are all those between the upper and lower broken lines, 
and additional storage is needed for the lower and upper triangular matrices too. 

2.5. Relative efficiency of the algorithms 

Examples of elapsed CPU time for the three methods are given in Fig.|4] Three NLTE cases of "physical interest" 
are illustrated here: carbon at 10 eV and 10'^ e.cm"^ with the 7 charge states and about 1800 levels, krypton at 1 
keV and lO''* e.cm"^ with Krxxiv-xxxiv ions and about 5600 levels included, and at 200 eV, 10^^ e.cm"^ with Krxix- 
xxvm ions and about 9700 levels included. In the last case, the kinetic system was too large to allow ourselves for 
a resolution by the Gauss algorithm in a reasonable time. It appears that while computation time increases rapidly 
with the matrix size, the CG and LU methods are orders of magnitude faster then Gauss elimination. Because as seen 
above the tested CG algorithm lacks of robustness, the preferred LU method is used in all the computations presented 
below. 



3. Configuration-average procedures and validity criteria 
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Configuration average (CA) procedures have been discussed in previous papers lllSL Il7l 12111 . It is the simplest 
and more natural way to considerably alleviate the matrix inversion task involved in stationary NLTE cases, and, 
interestingly enough, it usually provides rather accurate results. 

3.1. Configuration average(s) 

Assuming the ionic level / (resp. j) belongs to configuration a (resp. fi) the average rate from or to y6 is 

Rap^—^ giRij, with 8a = gi (2) 

iea iea 

where gi is the /-level degeneracy, and ga is the configuration degeneracy. We can set an alternate definition by 
including the Boltzmann factor 



Rap = ^ gi exp{-Ei/(kBTe))Rij Si exp (-EiKkBTe)) , 



(3) 



which can be assumed more "accurate" because the levels inside a given configuration usually have neighboring en- 
ergies and are likely to fulfill thermal equiUbrium. Of course, corresponding definitions hold for the average energies, 
energy rms, etc. 



Another option that generalizes both previous ones consists in arbitrarily defining some "configuration tempera- 
ture", i.e., 

Rap = gi expi-Ei/ikBT.Wij 2 gi exp (-Ei/(kBT,)) , (4) 

iea ' iea 

where Tc can be freely chosen. The dependence of the various physical quantities on the choice of TdTe will qualify 
the CA validity, as will be illustrated in the next sections, in cases where CA holds or does not hold. 

3.2. Validity criteria 

These CA formula being defined, it is important to elaborate a criterion that can properly qualify the relevance of 
this procedure, of course avoiding to go back to the time-consuming detailed level analysis. 

3.2.1. The natural criterion based energy dispersion 

The most intuitive criterion for the validity of CA may be written as 

< ^ Na^Ea « ksTe (5) 

a 

where is the a-configuration population (with N„ = I) and AEa the energy dispersion 

\l/2 



-VgiiEi-Eaf 



(6) 



This criterion is fairly easy to implement, since it only requires to solve once a simple CA kinetic system. However, 
as stressed before 11711 and checked once more here (see, e.g., subsection 15.3b . this necessary condition is far from 
being sufficient, a plain reason being that it involves a strongly averaged quantity, which furthermore depends on level 
energies and not on transition rates that can vary significantly inside a given pair of configurations. 

3.2.2. A criterion based on LTE test and its applicability range 

A criterion based on the comparison between the Saha-Boltzmann solution and of the solution of a partial rate 
kinetic equation has been proposed in previous papers iflTi I21I . It consists in checking how well the C A of a system 
containing rates that obey detailed balance — collisional rates in iflTl . any rates in 1I21]] — reproduces the Saha- 
Boltzmann solution. In more detail, one first evaluates using ^ average collisional rates R'^^ where the rates R' 
include only the collisional excitation, the collisional ionization, and the reverse processes. Then one solves the CA 
kinetic system 

that gives the configuration populations Na- Last, one compares this solution with the Saha-Boltzmann populations of 
the configurations A'^^. Since the R' rates obey detailed balance, the difference - A^^^^ is a measure of the validity 
of CA. This criterion has proven to give good results in a carbon plasma at 3 eV or 10 eV iIitIi . 
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However we must point some drawbacks of this approach. First, the criterion in its simplest form relies on a 
comparison of average net charges, which is known to be rather insensitive to the detail of the populations, regarding 
excited levels for instance. Second and more seriously, this criterion checks the CA validity close to LTE, and, if the 
studied situation is significantly far from LTE, it cannot be predictive. Examples in this direction wiU be provided in 
next sections. This criterion appears as necessary but not sufficient either. 

3.2.3. An alternate criterion based on fictive configuration temperatures 

A last approach consists in solving the CA kinetic equations using the various averages [Eqs. (|2]i, ([3]), (|4|] and to 
compare the results. The comparison may concern several physical quantities such as < Z* >, the radiative bound- 
bound rate 

^bb = Anions Yj P' Z^-^' ~ ^j^^'j 
' j 

Pi being the /-level population, £, - Ej the transition energy, and A,y the radiative rate, or the radiative bound-free rate. 
This criterion requires to solve at least two CA systems, and preferably three, one at Tc infinite, one at Tc = T^, and 
one at a fictive low e.g., 0.1 Tp. Examples will be given in the next sections. 

4. The krypton case: validity and limitation of CA 

The analysis of charge distribution and radiative losses in krypton presents a definite interest for various reasons. 
Krypton is used in numerous plasmas devices, and one must mention the availability of radiative loss measurements 
by Fournier et al jl'Jil . 

The present study deals with temperatures at least equal to 500 eV, and electron densities up to 10^^ cm charge 
states from XIX to XXXVII are accounted for. For such large temperatures, the collisional rates are low enough to 
provide a strongly non-LTE situation. The configurations included in the computation are listed in Table [T] In order 
to include an important set of excited states required for a pertinent description, this list is much more complete than 
the one considered for test purposes in Sec. |2] For instance at 1 keV, lO'"* e/cnr', the present computation includes 
45 927 levels vs 5 587 previously. The energy and rate computations in krypton as well as in tungsten are performed 
using the HULLAC suite 

This accounts for a large number of autoionizing configurations while keeping the NLTE computation tractable 
even in the detailed case. Lower temperatures would require a more complete set of ions but are within the reach of 
present method. 

The average charge, the radiative bound-bound (bb) losses and the radiative bound-free (bf) losses are plotted in 
Fig. |5]and|6]for electronic densities equal to lO''* and 10'^ cm respectively. Detailed computation and configuration 
average results are both displayed. Once again, the agreement between both computations is good. This is particularly 
true at r = 1000 eV and T = 5000 eV where the average charge is close to, respectively, 26 and 34, which corresponds 
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to the Ne-like and He-like closed shell ions. Then it is well known that all computations tend to give agreeing 
predictions. 



Since Fournier et al performed measurements of the radiative cooling coefficient in coronal Kr plasmas ll22ll . it is 
instructive to compare our 10'"* cm"^ data to theirs. The radiative cooling coefficient is defined for coronal plasmas as 

L'°' =< Z* > (Pbb + Ph + P«)INl (9) 

where Ps is the Bremsstrahlung loss rate, which we can derive from the approximate expression (formula 4.24.4 in 



Pff(erg/s/cm^) ^ 5.35 x 10"^"* < Z* > N^iT,)^^^ with A^, in cm"^ T, in keV. (10) 

More accurately, we should substitute < Z*- > / < Z* > to < Z* > in the rate ( fTOl ) but this looks unnecessary since 
this is just an estimation. At 1 keV and 10'"^ e.cm"^, the ff term is f ff ^ 1 .4 x 10^ cgs, less than the bf one and quite 
small with respect to the bb contribution. 

From Fig. 6 of Ref.|22|]. we get a radiative coefficient L'°' = 6 x 10"^^ W.m^ at 1 keV, while the present de- 



termination is 7.8 x 10"^^, in very good agreement if we consider the difficulty of this measurement. At 500 eV, 
the agreement deteriorates (1.6 x 10"^' SI measured, vs. 5.9 x 10"^' computed), while at 2 keV one gets again fair 
agreement (6 x 10"^^ measured vs 8.6 x lO^^*^ obtained here). 

A clear discrepancy between detailed and configuration average is visible on the bound-free losses at 500 eV 
and 10'^ cm"^: 2.82 x lO''* erg/s/cm"' in the detailed model, 4.43 x 10''* erg/s/cm"' in CA. Nevertheless this bf-loss 
ratio does not exceed 1.57, and bound-free losses are two orders of magnitude lower than bound-bound losses, which 
means that such a difference will hardly be noticed experimentally, only the sum being measured. However, one may 
try to look closer at the origin of this discrepancy. A first explanation could come from the averaging procedure that 
can include or not the Boltzmann factor (Eqs. [3] |2]l, the Fig. |6] corresponding to the case including the Boltzmann 
factor. The inspection of Table |3] demonstrates that the averaging procedure is not at stake. For all the plotted data, 
both configuration-average quantities would be indistinguishable at the drawing accuracy. One simply notices that the 
agreement between both averages increases with T^, as expected. 

A more thorough analysis arises from the detailed analysis of the radiative losses contributions. The configurations 
contributing the most to bound-free losses according to the CA calculations are listed in Table |2] One may see that, 
e.g., the \s^2s^2p^3p3d configuration of Kr^"*^ dominates for the bf losses in the CA computation (total population 
0.0147, contribution to bf losses 2.48 x 10'"* erg/s/cm^), while the corresponding levels contribute for a population of 
0.0219 and bf losses of 5.13 x lO'^* erg/s/cm^ in the detailed computation. 

Conversely, in the detailed scheme, the level contributing the most to bf losses is the Kr^^+ ground level (8.46 x 
10'^ erg/s/cm^), followed by the Ki^^^ {Is^ls^lp^ySoipiiiidsiiJ = 3 level (4.68 x lO'^) and Kr^^^ Ne-like ground 
level (4.65 x lO'^). 

As a rule one notices that a detailed inspection (e.g., at the configuration level) of populations or radiative losses 
may reveal large discrepancies between models that will remain hidden when only global quantities are considered. 



To check how the CA vaUdity criteria defined in subsection 13 .2. 3l behaves. we have done two additional compu- 
tations of the Kr plasma properties in configuration average. One was performed using the "infinite configuration 
temperature" scheme as defined by formula while the last one used a finite temperature (|4|l chosen below the 
electronic temperature Te/T,, = 0.1. The results of the various CA compared to the detailed computation are shown 
in Table [3] It appears that the average (|4]l agrees usually within 2% or better with (|2]l, except on bf losses at 500 eV, 
where the diff'erence exceeds 5%. This discrepancy is an indication of the (stronger) dilJerence between the CA and 
detailed values. A computation with Tc/Te = 0.01 would increase this diff'erence. 



5. The tungsten case: an example of configuration-average vaUdity breakdown 

Tungsten exhibits some interesting properties that qualify it for being one of the constituents of the divertor in 



modern tokamaks ll25ll . It presents a high melting point, a low erosion rate, and a low tritium retention. However, 



in such a plasma one may expect high radiative losses because of the richness of its spectrum. It also allows to 



benchmark X-ray spectrum calculations ll26ll . That is why tungsten was in the list of cases submitted to the NLTE-5 
workshop lfl6i . 

5.7. Description of the computation 

The list of configurations used in the present computation is similar to the one used in krypton (Table [Ij. The 
computations performed here are restricted to - 20 keV and 30 keV, at a density of Ne - 10'^ cm""*. The density 
lO^"* cm"^ was also considered, but since the present model does not account for continuum lowering (or pressure 
ionization), these results are not really significant and could be exploited only if compared to alternate models. The 
net charge states included are 59-72 and 61-74 for Te = 20 keV and 30 keV respectively (at 10^"^ cm"^, 60-74 and 
61-74 charges were considered). 

5.2. Average charge and radiative losses analysis 

Results are summarized in Table |4] Considering the < Z* > value, we check that detailed and configuration- 
average computations provide similar figures: at 20 keV one gets 65.004 and 64.547 respectively, and at 30 keV, 
67.389 and 67.265, where the CA is performed using Boltzmann ponderation as defined in (O. The bound-free radia- 
tive losses, mostly sensitive to the ground-state population of each ions, are also quite similar in both computations. 
However, a huge discrepancy is observed on the bound-bound radiative losses, since at 30 000 eV the computed 
figures are 4.866 x 10^ and 6.496 x lO'^* erg/s/cm-' in the detailed and CA cases respectively. 

This surprising feature cannot be explained by a drastic variation of the charge distribution, rather similar in both 

cases. To trace the origin of this discrepancy, one must resort to a finer analysis of the contributions of each pair of 

configurations to the bound-bound rate. Doing this, it appears that the huge bound-bound losses in CA case originate 

mostly from the transition between the ls^2s^2p and the ls^2s^2p configurations of the Wlxx ion. While the former 

configuration exhibits a plain fine-structure splitting of the 2p doublet (of about 1500 eV, which much less than the 
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electronic temperature), the latter one has a more complex structure. It contains 8 levels, the 5 lowest ones being at 
energies close to one of the ls^2s^2p doublet component, while the 3 upper levels are approximately 1700 eV above 
the 2/?3/2 upper component of the doublet. The two configuration levels therefore strongly overlap. Considering the 
radiative deexcitation rates, a detailed inspection demonstrates that the 5 higher levels of the ls^2s^2p configuration 
exhibit the largest radiative rates (above 10'^ s ' ), while they are noticeably the less populated. This fact is illustrated 
by Fig. |2l where the eight level populations ; e ls^2s2p^ as obtained from the detailed computation are plotted in 
ordinate while the abscissas stand for the individual radiative loss factor Tjjeis^is^ip^ijiEi - Ej). The anticorrelation 
between populations and rates appear strikingly on this plot. An additional explanation for the inadequacy of the CA 
computation in this case comes from the total population of the \s^2s2p^ configuration. The CA computation gives 
^ 4 X 10"^, while the sum of the eight level populations in the detailed approach is only ^ 3 x 10"^, lowering thus the 
radiative bound-bound losses by 5 orders of magnitude. The discrepancy on this configuration population originates 
also from significant variations in the coUisional and radiative detailed rates between configurations, and it has not 
been analyzed in more detail. 

5.3. About various criteria assessing the validity of configuration average 

One must now consider whether this strong discrepancy was expected, since the CA average remain much less 
computer-time consuming. 

First the plain criterion on the average energy dispersion is pointless here: one computes for the energy rms inside 
each configuration ponderated by its population ^ < AE >= 625 eV, far below = 30 000 eV. Then, the criterion 
of subsection l3.2.2l based on the comparison of the Saha-Boltzmann solution with the CA solution of a modified rate 
equation is of no help either The temperature being very high (and the density very low), the modified rate equations 
give a plain fully ionized plasma, the Z* < 74 ion population being less than 10"^. This solution is in excellent 
agreement with Saha-Boltzmann equation, but provides no relevant information on the real plasma, far from thermal 
equilibrium. 

One must then resort to the fictive-configuration-temperature criterion devised in subsection l3.2.3l The use of this 
criterion is illustrated in Table ID Therefore, in addition to the detailed computation, we have performed the three 
configuration averages proposed in Sec. [3] The first average (|2]i is the g ponderation, i.e., the infinite temperature 
limit; the "natural" average Q involving Boltzmann factors corresponds to the previously discussed CA value. For 
the "fictive configuration temperature" case (|4]i one has taken Tc/Te = 0.1, in order to test Tc well above and below 

Te. 

Considering first the Z* average, comparing the CA data obtained with an infinite temperature (column 3) with the 

data with low Tc (column 5, the column 4 being in between), one checks that these two figures are close. If we look 

now at the b-b losses (columns 6-9), we see that (|2|l, or (O which is close, differ by almost a factor of 2 at 20 000 

eV, and less but still significantly at 30 000 eV. This discrepancy between averages (O and must be interpreted 

as a failure of the configuration average, i.e., is a clue of the even much stronger discrepancy between the detailed 
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result and the ((Sj average. Once again, the interesting fact in using Q and Q is that this requires only two quite 
fast configuration average computations to check whether a detailed analysis is required. Moreover, looking at Z*, 
bound-bound and bound-free losses in this table, one notices immediately that only the bound-bound figure in CA 
formalism is highly dubious. 

In order to validate further the Tc criterion, we also performed computations with Te/T^ - 0.01. It appears that, 
at 20 000 eV, one gets < Z* >= 64.362 with Tc - 200 eV, not too different from the Ti-oo limit, which strengthen our 
confidence in this CA average. Conversely, in the same conditions, the bound-bound losses are 2.2 x 10^ erg/s/cm^ 
to be compared to the Tc-infinite CA value of 1.3 x lO'^ erg/s/cm^: This huge change clearly indicates a failure 
of the CA approximation when dealing with bb losses. Once again, the modest variation of the bound-free losses 
(8.48 X 10^ erg/s/cm"' to 8.60 x 10^ erg/s/cm"' when Tc varies from Jg/lOO to oo) insures that such observable is 
reasonably computed within CA approximation. 

6. Conclusions 

Reliable computations of plasma properties significantly out of local thermal equilibrium involving complex atoms 
require to solve large linear systems, for which various algorithms have been tested. Without claiming exhaustivity, 
it appears that the LU decomposition method ofifers a good compromise regarding stability, robustness, accuracy and 
computation speed. Alternate methods using, e.g., conjugate gradients with a preconditioner are worth consideration 
but certainly of less direct usage. We have investigated a second direction to simplify the kinetic system solution, 
based on configuration average. In addition to the already known criteria, that provide a necessary but not sufficient 
condition, we have tested a criterion based on the use of a fictive configuration temperature. This criterion proves 
to be very efficient including when dealing with situations far off' thermal equilibrium. Moreover, even negative 
configuration temperatures may be used, allowing a vast range of configuration average checking. Since in some 
situations, as illustrated by the bound-bound power losses in tungsten at high T^, this average performs poorly it is 
essential to use an efficient validity criterion and when necessary to dispose of a fast and robust algorithm for the very 
large kinetic systems that must then be solved. 
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Tables 



Ion 


N, 


Ni 


List of configurations 


KrXIX 


70 


3781 


[ls-2s^2p^']3s^3phtl, 3s^3p*3dnl, 3s3p'^'nl, 3p^3dnl\ ls^2s^2p^3s^3p^nl, 
ls^2s2p''3s^3p^nl 


KrXX-XXII 
{N = 4, 3, 2) 


63 


a 


lls^2s^2p^]3s^3p'^nl, 3s^3i^-^3dnl, 3s3p^*^nl\ \s^2p-2p^3p-3i^nl, 

\s^2s2p^3s^3p^nl 


KrXXin 


80 


6006 


[l/2i22/]3/3pn/, 3/3rfn/, 3s3p^nl, 3s3p3dnl, 3p\l, 3p^3dnl; 
ls^2s^2p^3s^3p^nl, ls^2s2p^3s^3p^nl 


KrXXIV 


75 


2239 


[l/2/2p«]3/nl, 3s3pnl, 3s3dnl, 3p3dnl; ls^2s^2p^3s^3pnl, 
l^2s2p^3sHpnl 


KrXXV 


61 


1499 


[\s^2s^2p^]3snl, 3pnl, 4/4/', l.v22.v22/753/3/'3/"; l.v2.v22/3/3/'3/" 


KrXXVI 


67 


2129 


li22/2/n/, \s^2s^2p^nln'l', l/2i2/3/3/'; li2i22p63/3/' 


KrXXVII 


83 


5397 


ls^2s'^2p^nl, l/2i2/n/, l/2/2/n/n'/', l/2i2p53/3/', 1/2^*3/3/'; 
li2/2/3/ 


KrXXVIII 


82 


6758 


ls^2s^2p'^nl, \s^2s2p^nl, ls^2s^2p^nln'l'; ls2s^2p^ 


KrXXIX 


74 


4739 


\s^2s^2p^nl, ls^2s2p*nl, \s^2p^nl, \s^2s^2p^nln' I' \ \s2s^2p\ \s2s2p^ 


KrXXX-XXXI 
(N = 2, 1) 


94 


b 


\s^2s^2p'^nl, l/2i2p"+in/, l/2p'^+2„;^ ls^2s^2p'^-^nln'l'; \s2s^2p^*^, 
Is2s2p^*^, mp^"-*, ls2s^2p^*^3l, U2s2p^*'^3l, \s2p^*'^3l 


KrXXXII 


135 


3568 


1/2/n/, l/2i2pn/, l/2p2„/, ls^2snln'l'; ls2s^2p^, Is2s2p^, ls2p*, 

isis lp3i, isZsZp 31, i-SZp 31 


KrXXXin 


80 


751 


\s^2snl, ls'^2pnl, Is^nln'l'; U2s^2p, Is2s2p^, ls2p^ , \s2s^3l, U2s2p3l, 
ls2p^3l 


KrXXXrV 


23 


106 


l^nl; \s2^, U2iP-, U2131' 


KrXXXV 


18 


59 


\snl\ 2/2/' 


KrXXXVI 


15 


25 


nl 



" Nl = 6988, 8980, 7319 for W = 4,3, 2 respectively 
^ Ni = 3607, 1915 for N = 2,\ respectively 

Table 1: Configurations included in Kr computation. Unless otherwise specified, one has n < 5,1 < n - I. When n 

and n' are involved in the same configuration, one has 6 < n + «' < 8 with I < n - < n' - \. Nc 'is the number 
of nonrelativistic configurations, Ni the number of levels. In the HULLAC computation, the configurations are usually 
separated in two groups within which interaction is fully accounted for; these groups are separated by a semicolon in 
the list, the second group containing excitation of an inner electron {Is, 2s or 2p). 
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Z Configuration CA: population CA: bf losses DL: population DL: bf losses 

(erg/s/cm-') (erg/s/cm^) 



24 ls^2s^2p''3p3d 0.0147 2.48 X lO^'' 0.0219 5.13x10 
23 ls^2s^2p^3s^3p 0.318 6.02 xlO^^ 0.0931 1.68x10 

25 15^25^2/35 0.157 4.50 xlO^^ 0.305 8.46x10 

Total 1 4.43 xlO^'' 1 2.82x10 



Table 2: Configurations with the largest contribution to the radiative bound-free (bf) losses in Kr at 500 eV and lO'* cm"^, in the configuration 
average (CA) case and in detailed level (DL) computation. The total population of the configuration and its contribution to bf losses (erg/s/cm^) is 
given in each case. 
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Te (eV) 




<Z* 


> 




Bound-bound losses (10'^ 


erg/s/cm-') 


Bound-free losses (10''* erj 


i/s/cm-') 


(eV) 


detailed 




(13JI 


/TTIi 
dill 


UcLallcU 






(14J| 


(JcLdllcU 


(UP 




dlJl 


son 




7^ 7zlQ 






5.949 


9.327 


9.332 


Q ^81 

JO i 


2.823 


4.406 


4.431 




1000 


25.786 


25.654 


25.654 


25.650 


1.919 


2.333 


2.331 


2.318 


2.465 


2.495 


2.499 


2.532 


2000 


28.457 


27.826 


27.830 


27.860 


2.553 


3.049 


3.053 


3.097 


3.299 


3.095 


3.097 


3.110 


5000 


33.297 


33.234 


33.234 


33.234 


0.9903 


1.347 


1.345 


1.319 


4.638 


4.573 


4.573 


4.573 



Table 3: Kiypton plasma properties at lO'** e.cm The average net charge < Z* >, radiative bound-bound and bound-free losses are given in the detailed-level models (columns 2, 6, 10) and in 
the configuration average approximation. This average is performed using in the Boltzmann factor a configuration temperature Tc infinite (Eq.|2j, equal to the electronic temperature Tg (Eq.[3) or 
such as Tc = O.ir^ (Eq.gl. 



Te (eV) 
(eV) 


<Z* > 

detailed © Q (Hi 


Bound-bound losses (10^ erg/s/cm"') 
detailed © © © 


Bound-free losses (10^ erg/s/cm^) 
detailed © © © 


2x lO'* 
3 X 10" 


65.004 64.541 64.547 64.571 
67.389 67.263 67.265 67.254 


6.370 1.290[6] 1.223[6] 6.766[5] 
4.866 6.683[6] 6.496[6] 4.866[6] 


9.238 8.601 8.606 8.635 
9.917 9.396 9.399 9.424 



Table 4: Same as table[3]but for a tungsten plasma at lO'* e.cm ^ . The notation 1.290[6] stands for 1.290 X 10*. 
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Figure 1: CoUisional-radiative rate matrix structure: an example in krypton at = 1000 eV and N,, = lO'* e/cm'. The number of equations 
is here 5586. The net ion charge is indicated for each diagonal bloc. The cumulated number of levels is indicated in abscissa and ordinate. The 
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Figure 2: Average charge (upper figure) and its variation from one step to tlie previous (lower figure) in Kr as computed witli conjugate gradient 
method stalling on Saha-Boltzmann populations. The horizontal red line coiresponds to the average charge as determined by the stable Gauss 
algorithm. The electron temperature and density are 1000 eV and lO''' cm"^ respectively. The e parameter of Eq. Q is 10"". 
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Figure 3: Same as Fig.[2]but starting on zero populations. 
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Figure 5: Average charge, radiative bound-bound and bound-free losses for a krypton plasma at lO'^ cm ^ electronic density as a function of the 
temperature in eV. Computations are peiformed using detailed levels and configuration average. 
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Figure 7: Radiative losses from and populations of the eight levels belonging to the l.v^2.v2p~ configuration of W^'^*. The lower levels j of the 
radiative decay belong to the l.v^2s^2/) configuration. The detailed level computation is done at 30 000 eV and lO'"* e/cm'. 
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